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We develop ion-ion pair potentials for Al, Na and K for densities and temperatures relevant to the warm-dense- 
matter (WDM) regime. Furthermore, we emphasize non-equilibrium states where the ion temperature Ti differs 
from the electron temperature Te. This work focuses mainly on ultra-fast laser-metal interactions where the 
energy of the laser is almost exclusively transferred to the electron suh-system over femtosecond time scales. 
This results in a two-temperature system with Te > Ti and with the ions still at the initial room temperature 
Ti = Tr. First-principles calculations, such as density functional theory (DFT) or quantum Monte Carlo, are as 
yet not fully feasible for WDM conditions due to lack of finite-T features, e.g. pseudopotentials, and extensive 
CPU time requirements. Simpler methods are needed to study these highly complex systems. We propose to 
use two-temperature pair potentials Uuir, T, Te) constructed from linear-response theory using the non-linear 
electron density n(r) obtained from finite-T DFT with a single ion immersed in the appropriate electron fluid. 
We compute equilibrium phonon spectra at Tr which are found to be in very good agreement with experiments. 
This gives credibility to our non-equilibrium phonon dispersion relations which are important in determining 
thermophysical properties, stability, energy-relaxation mechanisms and transport coefficients. 
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1 Introduction 

Warm dense matter (WDM) is a region of the phase diagram located at the conjuncture of the plasma state and 
the condensed-matter state. The density is high enough that ions are strongly correlated and can be found in a 
mixture of different ionization states. The ion temperature Ti can be different from the electron temperature Te 
which can be as high as the Fermi temperature or even more, and hence the electrons can be partially degenerate. 
Such WDM states are of great interest at the fundamental level, but also for their technological applications 
such as laser ablation IT], inertial-confinement fusion 121, Coulomb explosion El, etc., and are also of relevance 
to astrophysics Elia. It is a great challenge to predict the two-temperature quasi-thermodynamic properties 
or transport and radiative processes of such basically non-thermodynamic states. The interest for WDM in the 
theoretical physics community has greatly increased since the advent of the possibility of creating these extreme 
conditions in the laboratory. However, the experimental techniques required for this purpose have the crucial 
particularity of always producing two-temperature WDM with Ti ^ Te, which is problematic in the current 
theoretical framework. The present work focuses on the case where Te > Ti, which can be observed when a 
femtosecond laser pulse interacts with a thin metal foil 017]. Since the electron-ion mass ratio is very small, 
the laser energy is almost totally transferred to the electron system, thus increasing the electron temperature 
from the equilibrium initial temperature to Tg. Meanwhile, since the transfer of energy from the laser to the 
electron subsystem, as well as the equilibration of the latter, happens much faster than the electron-ion energy 
relaxation, the ion temperature Ti remains at room temperature Tr with the ionic-lattice structure staying intact as 
the ions have no time to move. Using a two-temperature model (TTM) 0|9] to describe the electron-ion energy 
relaxation, it has been determined that the electron-ion relaxation time Tei is of the order of picoseconds, while 
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Tee and Tii are orders of magnitude shorter. The challenge then is to develop theoretical tools designed for such 
non-equilibrium and strongly-correlated systems. 

Within the framework of density functional theory (DFT) ifTOl . the total free energy of a system of electrons 
and ions T'[n(r), /o(r)] can be written as a functional of the one-electron density n(r) and the one-ion density 
p(r). The energy minimization results in a set of two coupled equations for n(r) and p(r). The first set is the 
integral equations for the ion distribution p(r) which obeys classical physics, and hence can be solved using the 
hypernetted-chain (HNC) approximation within the Ornstein-Zernike (OZ) equation if the ion-ion pair potential 
Uii (r) is known. The second set is the Kohn-Sham HD equations for electrons in the external potential caused by 
p(r). The HNC approximation, or its modified version with bridge corrections, provides the correlation potentials 
for the ion subsystem which does not have exchange effects. However, a key element in the DFT-equations of 
the electron subsystem is the exchange and correlation (xc) free energy functional and its density 

derivative which is the xc-potential. Major DFT codes — e.g. ABINIT ifT^ — are designed for ground-state 
calculations using ground-state adapted pseudopotentials and do not include many finite-T features. They offer a 
large variety of ground state xc-functionals Exc\n] but do not currently implement finite-T xc-functionals. Since 
DFT calculations can be quite sensitive to the choice of Exc[n] in certain regimes, it seems very important to use 
a finite-T functional Fxc\n^ Te] to adequately study WDM systems. Furthermore, in the DFT extension to finite- 
T by Mermin CD, it is required to include the Fermi distribution f{Ei,Te) into the finite-T electron density 
n{r,Te) calculation. Since Tg can be very high in the WDM regime, it requires many excited Kohn-Sham- 
Mermin eigenfunctions in each part of the self-consistent procedure which considerably increases the calculation 
time. Moreover, standard DFT simulations are designed for crystal or molecular structures, which excessively 
limit the range of T^ accessible to study WDM. 

Higher ion temperatures can be studied using molecular dynamics (MD) simulations. MD is a powerful 
method to study the ion system in the liquid, gaseous or plasma phases. It has been highly parallelized to make it 
possible to simulate very large systems. The main input of MD simulations is the ion-ion pair potential which is 
required for each interaction in the system. However, most pair potentials, as implemented in comprehensive MD 
packages such as LAMMPS m, are semi-empirical constructions. A particular example is the embedded-atom 
method (EAM) lfT9ll . fitted to reproduce particular properties under ambient conditions while failing at predicting 
other ones, especially those of interest in two-temperature WDM. 

In the following, we propose to construct ion-ion two-temperature pair potentials (TTP) Uu (r, Ti,Te) based on 
fundamental considerations. The first part of this paper explains how we construct the non-equilibrium electron- 
ion pseudopotentials Uei{r, rs,Tei), which depend on both temperature Tei = {Te, T) and density via the elec¬ 
tron Wigner-Seitz radius Vg. The most important quantity to compute is the free electron density nf{r,Tei) 
around an ion of effective charge Z. This density is then used to construct the pseudopotential from linear- 
response theory, going beyond the random phase approximation (RPA) by adding a finite-T local-field correction 
(LFC). The second part explains how to construct the pair potential Ua from Uei and shows the importance of 
the LFC. The total ion-ion interaction is obtained by adding the direct-ion Coulomb interaction to the indirect 
interaction through the electron subsystem. In the last part, we compute the equilibrium and non-equilibrium 
phonon spectra uj^{k, Tei) of three simple metals to test our TTP Uii{r, Tei). We examine aluminum as well as 
sodium and potassium which are typical free-electron metals under the conditions studied. 


2 Non-equilibrium electron-ion pseudopotentials 

According to the original formulation of DFT by Hohenberg and Kohn ifTOll . the electron density n(r) rather than 
the many-body wavefunction is considered to be the fundamental quantity describing a quantum system so that 
all physical properties of interest can be expressed as a functional of n(r). The self-consistent scheme proposed 
by Kohn and Sham HD enables one to reduce the many-body system to an effective “one-electron” problem by 
putting all the missing information about the electron-electron interactions in the xc-functional Exc[n]. However, 
since electrons in the WDM regime are far from zero temperature, it is crucial to use the extension of DFT to finite 
temperatures. Furthermore, one could wonder if using the DFT for systems with Ti ^ Te is in contravention with 
the fundamentals of DFT, but since calculations are done for a single fixed ion surrounded by its T^-dependent 
density which is immersed in an electron liquid at Tg, the applicability of DFT is satisfied. 
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2.1 Electron density from finite-T DFT 

In extending the work of Hohenberg and Kohn to finite-T, Mermin showed that the electron density defines 
uniquely the external potential V(r) of the system lfT3l . The Mermin-Kohn-Sham equations reduce the finite-T 
many-electron problem to a one-electron problem, giving the interacting electron density n(r, Tg). 

^ 1 
n(r,T) = nb(r)-fn/(r) = ^/(ei,T)|</)i(r)|2, where /(e»,T)= ^ _ (1) 

with nt,(r) and n/(r) the bound and free electron density, f(ei) the Fermi occupation function, ijP{T) the non¬ 
interacting chemical potential, /3 = 1/ k^T the inverse temperature (with /ch = 1 in our units), and (pi (r) the 
Mermin-Kohn-Sham wavefunction of eigenvalue e^. At T = 0, /i° is the Fermi energy ep- As mentioned in 
the introduction, the presence of the Fermi occupation factors, and the need for high angular-momentum states, 
greatly increase the CPU time since the convergence of the energy is quadratic in the number of bands needed. 
One way to decrease the simulation time and make it accessible to WDM conditions is to reduce the number of 
electrons in the system by introducing the effective ion charge Z which appears in the pseudopotential. However, 
the formulation of finite-T pseudopotentials remains largely unexplored, and standard codes simply use the zero- 
T ones which are inadequate for many WDM applications Ha. 

2.1.1 Effective charge Z and the neutral pseudoatom approximation 

The main objective is to replace the nuclear charge Z by an effective charge Z = Z — rib where rib is the 
number of bound electrons. This is calculated by integrating the bound electron density nb(y) over the ion sphere 
iniiisi. The main difficulty is to determine how to split the total charge density into a bound and a free part 
n(r) = ribiv) -f n/(r) since bound electrons may extend beyond the Wigner-Seitz radius r^,. In this study, 
there is no ambiguity in the value of Z for Al, Na and K in the given density and temperature range. Using Z, 
one may attempt to construct an electron-ion pseudopotential which reproduces the free electron wavefunctions 
outside a core radius r^- However, a large number of continuum wavefunctions are needed in the WDM regime 
associated with finite-T problems na, and a large number of ions have to be included in the simulation box of 
the usual type of DFT calculations, increasing significantly the CPU time. We simplify the problem using two 
key ideas from DFT. (i) Only one nucleus at the origin is considered, and spherical symmetry can be exploited; 
hence r will be replaced by r when appropriate. The charge density n(r) around the central nucleus, rather than 
its wavefunctions are considered, (ii) Instead of including every ion at their locations we consider only the 
one-ion charge distribution p(r) = J^i ~ u)- Here we follow the correlation-sphere-DFT model of Dharma- 
wardana and Perrot and compute the electron density around only one ion of effective charge Z surrounded by 
the self-consistently calculated spherically averaged ion density distribution ufrl llSOl . Once the free electron 
charge density n/(r) at an ion in the given WDM environment is determined, a simple local pseudopotential can 
be extracted as explained in sec. 12.21 

For the elements that we study in this work, it is possible to simplify the calculation of n/(r) even more 
by using the neutral-pseudoatom (NPA) approximation El. The NPA approach is valid for solids, liquids and 
plasmas at low or high densities and temperatures as long as Z can be clearly defined so that core-electrons 
are not important in the ion-ion interaction. Thus, the ion density is replaced by a uniform positive background 
(jellium) except for a spherical neutral cavity of radius rws (ion Wigner-Seitz radius) around the central ion. In 
this model, the ion temperature Ti is merely a parameter which determines the value of rws from the density of 
the material at Ti. All these simplify greatly the solution of the Kohn-Sham-Mermin equations. 

2.1.2 The exchange and correlation free energy functional 

To compute n(r) around the single ion, the finite-T xc-potential given by Perrot and Dharma-wardana is used 
in the Kohn-Sham-Mermin equations. It has been constructed as follows. The fully non-local xc free energy 
per electron f^c = F^c/N is given by the integral of the spin-dependent electron-electron correlation function 
hijir) = py (r)-l over the coupling constant A 

^ f ^ A) -f %(r, A)]. (2) 
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Using the classical-map hypernetted-chain (CHNC) ED method, hij (r) and fxc have been calculated by Perrot 
and Dharma-wardana for different electron densities Tg and temperatures Tg. Their results appear to be in good 
agreement with the T = 0 and finite-T quantum Monte Carlo (QMC) calculations E2ll . thereby supporting 
the use of the CHNC approximation for problems where QMC data are not available. The method consists in 
mapping the quantum electron-fluid onto a classical Coulomb fluid using an effective electron-electron potential 
and an effective electron temperature to simulate quantum effects. It makes it possible to solve the OZ integral 
equations with the HNC approximation so as to get the finite—T interacting (r, A, Tg). The fxc{fs,Tf.) 
obtained this way has been parametrized to a simple form ll23l which makes it very practical. 

2.2 Electron-ion pseudopotentials from linear response theory 

Once the effective charge Z and the free electron density n / (r) have been calculated for a particular and Tei, it 
is possible to construct a weak local (s—wave) electron-ion pseudopotential Uei{r) using linear response theory 
ED- In g-space, it takes the form 


Ueiisii ; ^ei) — j'^si ^e): 


(3) 


with Xee the finite-T interacting electron linear response. Only the q 3kp region is relevant as the ion-core 
is excluded. To go beyond the RPA, we include a local-field correction Gq so that Xee is given by 


Xeeiq,rs,Te) 


_ xo{q) _ 

l-Vq{l-Gq)xo{q) 


with 



G, 



(4) 


Here Xo{q) is the non-interacting finite—T Lindhard function which may use an effective electron mass m* if 
needed E4ll . The EEC is used in the local-density approximation (EDA) and is evaluated from the ratio of the 
non-interacting and interacting finite-T compressibilities 70 and 7 respectively of the electron subsystem. The 
interacting compressibility can be evaluated from the finite-T xc-potential. It is important to note that the NPA 
approach is valid if resulting pseudopotentials satisfy Uei{q) /{— ZVq) < 1. 


3 Two-temperature ion-ion pair potentials 

3.1 The ion-ion pair potentials 

It is now possible to construct two-temperature pair potentials by adding the direct ion-ion interaction to the 
indirect interaction through the electron fluid such that 

Uii{q,rs,Tei) = Z‘^{Ti)Vq + \Uei{q,rs,Te^)\'^Xee{q,rs,Te). (5) 

The extension to a mixture of ions of different types or different ionizations is discussed in Ref. Il25l . In the 
present work, we present two-temperature pair potentials resulting from this procedure for three different simple 
metals. We examine aluminum, which is the “standard” free electrons metal with three valence electrons {Z = 3), 
followed by sodium and potassium, which are alkali metals with one valence electron (Z = 1). Since most laser 
pulse-probe experiments are done at room temperature, we computed pair potentials with Ti = = 0.026 eV 

whereas Tg is increased up to some value, e.g. 6 eV, that is experimentally relevant. The resulting TTP for Al, Na 
and K are presented in Eig. [1] The pseudopotentials used to construct these TTP are fitted to a Heine-Abarankov 
form (HA). A simple discussion of HA may be found in Shaw and Harrison E6l . The local HA pseudopotentials 
use a core radius Vc and a constant core-potential A for r < r^, while the potential for r > rg is Z /r. The 
equilibrium lattice parameter at room temperature ap resulting from such potentials is usually 1-3% in error. We 
have slightly modified the rg and A parameters to obtain nearly the correct ap. This is justified since the ion 
positions are held fixed at their room-temperature values in the WDM studied here. 

An important characteristic of these pair potentials is the presence of long-range oscillations which result from 
the Eriedel oscillations in the electron density. These are related to the discontinuity in the Eermi distributions 
in WDM and zero-T systems. They are clear non-linear effects included in our model from the self-consistently 
calculated electron density and exist even in linear response theory which is used to construct the pseudopoten¬ 
tials. Since even the first minimum is basically a Eriedel oscillation, it is very important to accurately calculate 
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Fig. 1 Two-temperature pair potentials for Al, Na and K. Ions are at T; = 0.026 eV and the electrons at Te = 
0.026 eV (black), 1.25 eV (purple), 2.00 eV (blue), 3.00 eV (red), 4.50 eV (orange) and 6.00 eV (gray). 


their amplitude. It has been shown that fitting pair potentials to elastic constants or phonon dispersion curves 
results in non-unique amplitudes lIZTl . We can see that the screening effect depends greatly on the number of va¬ 
lence electrons as the amplitudes of oscillations in the Na and K potentials damp much faster than in the Al case. 
Furthermore, as increases, the electron screening becomes less and less Friedel-like and the pair potentials 
exhibit only the repulsive behavior typical of Yukawa potentials. At such Tg, the crystal structure usually looses 
its stability, resulting in a Coulomb explosion or melting of the metal, depending on the heating time scales. 


3.2 The importance of the local-field correction 

To illustrate the importance of going beyond the RPA form of the screening function x, we computed the pair 
potential for Al at Tei = with and without the LFC; the results are presented in FiglJ] 



Fig. 2 Comparison between Al pair potential at Tei = Tr obtained using the RPA (black line) and with the 
LFC included (blue line) in the screening function. 

We observe a notable difference between the two pair potentials. The RPA-potential does not exhibit the first 
minimum shown by the LFC potential. This difference is crucial since this first minimum defines the nearest- 
neighbour position in the crystal; thus the RPA does not give a realistic description of solid Al. 


4 Phonon spectrum of laser-heated simple metals 


Using the TTP, it is possible to compute the phonon energies by constructing the dynamical matrix Il28ll : 


D(q) = ^D(r,)e-*‘^--- 

i 


with 


_ d'^Uujri - 
du{ri) clu(O) 


( 6 ) 


Here u is the small displacement in the position of ions at and at the origin 0. Phonon frequencies are then 
given by a;^(q) = \/A^(q)/M, with A^(q) the eigenvalues of the dynamical matrix. Using pair potentials, it 
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is extremely rapid to compute the phonon dispersions for any crystal structure. To validate our simple model, 
we compare it with a more fundamental linear-response method using the ABINIT package ifTSIl . The latter is, 
strictly speaking, only applicable to systems in equilibrium, with Ti = T^. Determining phonon spectra in DFT 
is usually done using the response function ll29l [30l method where the first derivative of every wavefunction is 
computed. This requires even more CPU time. Typically, while the ABINIT phonon calculation for aluminum 
may take many hours on a computer cluster, the pair potential approach is nearly ‘instantaneous’ on a laptop. 

The ions in the ultra-fast laser-heated solid are not at their equilibrium positions but frozen in their ‘old’ 
positions. The phonon description assumes that the elementary excitations are oscillations of the ions about their 
equilibrium positions, with the first derivative of the potential, i.e., the force, equal to zero. This is increasingly 
invalid as exceeds Ti. However, this ‘harmonic approximation’ is assumed in both the ABINIT and the pair- 
potential phonon calculation given here. In reality, the first derivative (i.e., the Hellman-Feynman force on each 
ion) is large and will lead to an explosive acceleration of the ions outwards rather than generating harmonic 
oscillations. A consideration of the time scale for one oscillation, and the time scale for the motion of an ion by 
about a wavelength, suggest that only the higher frequency modes can have a physical meaning. Nevertheless, we 
present the full phonon spectra in order to compare with similar efforts of other workers (using ab initio codes) 
who have provisionally assumed the existence of these non-equilibrium phonons. Besides the first-derivative 
being negligible, the second derivative of the total potential, being the elastic constant for the harmonic motion, 
must be positive at the ion positions. The higher derivatives give anharmonic effects. Thus, at sufficiently high 
Te, the phonon frequencies become purely imaginary, as the pair-potential becomes more like a Yukawa potential. 
It is extremely important to experimentally determine ‘phonon spectra’ of WDM systems to establish the scope 
and validity of the phonon description of the elementary excitations of these two-temperature solids. 


4.1 Aluminum 

At room temperature, aluminum exhibits a FCC crystal structure with an experimentally determined lattice pa¬ 
rameter ttAi = 4.05A which is used as an input (via the density) into our TTP construction. On the other hand, our 
ab initio simulations are done using the LDA with the norm-conserving pseudopotential method. We used a 165 
eV kinetic energy cutoff and a 20 x 20 x 20 Monkhorst-Pack k-point grid to compute the planewave expansion of 
the pseudowavefunctions. The ground-state calculation gives us an equilibrium lattice parameter = 4.0lA 
and phonon frequencies, at T^i = = 0.026 eV, in very good agreement with our TTP model and with experi¬ 

ment ED, as presented in Fig. 3. We computed the phonon spectra for the laser-heated situation where ions are 
kept ‘frozen’ at their initial equilibrium position, and with Ti = T^, whereas the electron temperature is increased 
up to Te = 6 eV. The non-equilibrium phonon spectra obtained with our TTP model and with ABINIT-DFT, by 
simulating more electronic band and using the Fermi-Dirac occupation function to mimic finite-T effects , are 
also displayed in Fig. 3. 
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Fig. 3 Equilibrium (full lines) and non-equilibrium 
(dashed lines) aluminum phonon dispersion relations 
computed with the TTP model (blue) and ABINIT 
(black), compared with experimental data (circles). 


Fig. 4 Phonon frequencies of laser-heated aluminum 
at Ti = Tr A Te for longitudinal (circles) and trans¬ 
verse (squares and diamonds) modes at symmetry points 
X (black), L (blue) and K (red). 
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As Te increases, we find that both methods show phonon hardening. The TTP model predicts a maximum 
frequency increase of 29% at the symmetry point X passing from 9.54 THz to 12.32 THz for the longitudinal 
mode. As mentioned by Recoules et al 1^ . the transverse mode frequency at the symmetry point L is not really 
affected, passing from 4.43 THz to 4.62 THz. Nevertheless, our results clearly illustrate that all other branches 
are also modified by the laser-heating mechanism. At room temperature, the TTP model predicts frequencies 
lower than DFT calculations whereas at Te = 6 eV they are higher. As Tg increase, finite—T effects become 
more and more important; further studies are needed to understand if this difference comes from the finite—T 
features which are absent in the ABINIT-DFT simulations. The fast computation of TTP by this method permits 
us to compute phonons for many different Tg as presented in Fig.4. The ability to calculate phonon spectra ‘on the 
fly’ can be of great utility in determining electron-ion energy relaxation and non-equilibrium transport properties 
such as the thermal and electrical conductivity. 


4.2 Sodium and potassium 

Sodium and potassium are alkali metals having a BCC crystal structure at room temperature with experimentally 
determined lattice parameters of UNa = 4.23A and ok = 5.23A, respectively, which will be used as input in the 
TTP construction. The ABINIT-DFT simulations have been done within the generalized-gradient approximation 
(GGA), a 24 X 24 X 24 Monkhorst-Pack k-point grid and a 1360 eV kinetic energy cutoff for both Na and K. 
Ground state calculations result in lattice parameter = 4.20A and = 5.25A, respectively, and yield 
phonon frequencies in agreement with our TTP model and with experiment ll^l34l as presented at Fig.5. 



Fig. 5 Equilibrium Tg = Ti = Tr (full lines) and non-equilibrium Te = 6 eV (dashed lines) phonon dispersions for a) 
sodium and b) potassium computed with the TTP model (blue) and ABINIT (black), compared with experiment (circles). 


Once again, our TTP results are in excellent accord with experiment for the equilibrium case which gives us 
confidence in applying the method to the non-equilibrium situation. In the case of Na and K, we observe that most 
frequencies increase but also that the shape of the dispersion relation changes as Tg reaches 6 eV. We can see new 
maxima appearing between symmetry points T — H and between 7T — P. In both cases, the maximum increase 
occurs at the point N with an augmentation for the frequency of the longitudinal branch of 37%. Finite—T ab 
initio simulations for these metals are still not completed. Hence the comparison of the predictions of our TTP 
model with more fundamental calculations is not yet final for this case. 


5 Conclusion 

As WDM studies are becoming more and more topical and important in many fields of physics, we have addressed 
the problem of describing theoretically such highly-correlated systems via a simple model based on transparent 
physics. Furthermore, to compare predictions with experiment, it is necessary to describe non-equilibrium sys¬ 
tems where Ti ^ Tg. Since ab initio simulations can be very time consuming in this regime, we have developed 
simple two-temperature pair potentials as the ‘work-horse’ of the method. They are constructed from fundamen¬ 
tal considerations by computing the electron density from finite-T DFT and then using linear response theory 
with a finite-T local-field correction. Pair potentials resulting from this procedure produce long-range interac¬ 
tions resulting from the presence of Friedel oscillations in the electron density. These TTP were used to compute 
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phonon dispersion relations for equilibrium Al, Na and K. The excellent accord with experiment and DFT sim¬ 
ulations (using the ABINIT code) establishes support for the validity of our approach for the non-equilibrium 
case. In all cases, our model predicts an increase of the frequencies, which is different for each mode branch, 
supporting the phonon hardening hypothesis. At higher Tg, where thermal effects are more important, frequen¬ 
cies obtained for Al with the TTP model are found to be higher than the ones obtained by DFT simulations. This 
discrepancy might come from the non-implementation of a finite-T exchange and correlation energy functional 
and the lack of other finite-T features in ab initio codes like the ABINIT. To validate this hypothesis, further work 
will be done to implement our finite—T xc-functional in a major DFT code. The possibility to quickly compute 
TTP for many different combinations of T^ and Tg could also be used to predict electron-ion energy relaxation 
and non-equilibrium transport properties. Another important application of our TTP is to incorporate it in laser 
ablation simulations where pair potentials dependent on both temperatures are needed. Finally, construction of 
two-temperature potentials for more complex metals, such as gold and copper, would be of great interest since 
many WDM experiments llTSl have been carried out with these metals. 
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